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r1 Abstract 

2 Nonlinear ODEs invariant under the group SL(2,R) are solved numerically. We show that 

^ solution methods incorporating the Lie point symmetries provide better results than standard 

'~~' methods. 
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1 Introduction 

Historically Lie group theory started out as a theory of transformations of solutions of ordinary 
and partial differential equations. Differential equations are still one of the most important appli- 
cations of Lie groups [6|. The most common use of Lie group theory in this field is to perform 
symmetry reduction. For ordinary differential equations (ODE), this means that Lie point sym- 
metries are used to reduce the order of the equation. If the symmetry group is large enough the 
order of the ODE can be reduced to zero. This is equivalent to obtaining the general solution of 
the ODE, possibly in implicit form. For partial differential equations (PDF), symmetry reduction 
means a reduction of the number of independent variables and usually leads to exact analytical 
solutions, albeit particular ones. 

The purpose of this article is to discuss a different application of Lie groups in the theory of 
differential equations and to present some new results and new examples. This application can be 
called symmetry preserving discretization of differential equations and its purpose is to improve 
numerical methods for solving differential equations. In this article we restrict ourselves to ODEs. 
For recent reviews with references to original articles see L5,2]- 

The idea is to start from a given ODE of order N 

F{x,y,y',y",...,/'^^) = (1.1) 

and its known Lie point symmetry group G with a Lie algebra L (the symmetry algebra) reaUzed 
by vector fields of the form 

X = ax,y)d, + cl>{x,y)dy (1.2) 

( |1.1| ) is replaced by an invariant difference scheme, i.e. a system of two equations 

Ea(n, Xn+K, Xn+K+l, ■■; Xn+L,yn+K,yn+K+\, —,yn+L) = (1.3) 

a = l,2 L-K^N, N>M 
relating the variables x and y in M different points. 

The scheme ( L3 1 is constructed so as to be invariant under the same group G as the ODE ( |1.1[ ). 



This means that equations ( 1.3 1 must be annihilated on their solution set by the prolongations of 



the vector fields (1.2 1 



prXiEa)\E,=E2=0 =0 (1.4) 

where the discrete prolongation is 



n+L 



prX= 2 {^(Xi,yi)d,, + 4>(xi,yi)dy,] (1.5) 



i=n+K 

In practise, this means that we can proceed as follows: 
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1. Find the N-th order differential invariants I%x,y,y', ..^y'-N)), j = 1, ..., J, of the group G and 
rewrite the ODE ( |1.1| ) in terms of these invariants : 

Fir,,q,...,I'j) = (1.6) 

2. Find the difference invariants Ia{n, x^+k, ■■■, Xn+L,yn+K, ■■■,yn+L) of the same group G. 

3. Expand the difference invariants in Taylor series about some reference point, say (x„, j„). 
Choose such difference invariants that we have 

Ij = I'j + 0(e) (1.7) 

i.e. such that the leading terms in the expansion coincides with a corresponding differential invari- 
ant (6 ^ is the continuous limit). 



4. Write the difference scheme ( |1.3| ) in terms of difference invariants Ij. In the continuous limit 
we will have 

Ei(/[,...,/5) = ^° F(/[,/2^...,/5) = (1.8) 

E2(/;^,...,/^) = ^ = 

In the invariant discretization the lattice is not given a priori, but emerges as part of the solution 
of the difference scheme. 

In physics and other fields of science, symmetries of a system are often better known than the 
dynamics and symmetries of equations can be more important than explicit solutions. Preserving 
symmetries in numerical calculations can be expected to improve the results, specially global fea- 
tures of solutions. 

The invariant discretization should be compared to what we will call standard discretizations. 
The lattice is given a priori, usually a uniform one. The ODE ( |1.1[ ) is discretized by replacing all 
derivatives by usual finite differences, e.g., 

m(x„+i) - m(x„) 
Ux = , Xn-nn + xo (1-") 

and similarly for higher order derivatives. 

Recent articles devoted to theoretical aspects of invariant discretization of ODEs include EUll 
m. In Ref.Ill it was shown on several examples that the symmetry preserving schemes provide 
better accuracy than standard ones and numerical solutions close to singularities where standard 
schemes fail. Here we will show that similar results hold in other situations. 

2 A two dimensional realization of the algebra sl(2,R) 

2.1 The Lie algebra and the invariant ODEs 

Let us consider the sl(2,]R) algebra with a basis realized by the vector fields 

Xi = dy, X2 = xdx +ydy, X-i = Ixydx +y^dy . (2.1) 
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It can be extended to a gl(2,]R) algebra by adding 

X4 - ydy 



(2.2) 



to the basis. Let us prolong these vector fields so they act on functions F(x,y,y',y",y'"). The ac- 
tion of the corresponding SL(2,R) group on the prolonged space with local coordinates {x, y, y',y",y' 
allows two differential invariants, namely 



^ 2xy" + y' 



n = 



,, x^iy'y'" - 3y"^) 



u'3 



r - 
h - 



u'5 



7- - y- 

Using these invariants we can write a second and a third order invariant ODE, namely 



/3 



and 



2xy + y = yy 



x\yy"-3y"^) 

— y-5 = ^(A) 



(2.3) 



(2.4) 



(2.5) 



where 7 is a constant and F(z) is an arbitrary function. The ODE ( |2.5| ) will be invariant under the 
group GL(2,R), including the dilatations generated by X4 if we restrict F{z) to be F{z) - az^^^- 
Eq. (|2.5|) specializes to 



2/ / /// 



.."2^ 



x\y'y"' - 3/'") = a{2xy" +y'fl^y'^l^ 
where or is a constant. 



(2.6) 



By construction the ODEs ( |2.4| ) and ( |2.6[ ) have symmetry algebras that make it possible to re- 
duce them to quadratures. 



For eq. ( |2.4[ ) this provides two explicit solutions 

2 



y\2i^) = ' 



yh±- ^C-yx CtO 

yh + —px C = 0. 

Vr 



(2.7) 



The two branches of the solution for C t^ intersect for x = C/y where we have y\ = yj = yt- 
After they become complex, the solution for x = C/y remains finite but all its derivatives become 
infinite. 



For eq. ( |2.6[ ) the quadratures lead to an implicit solution : 



y=yo + C[ I e'o 
Jo 

where f(x) satisfies 



fis)ds 



dt 



(2.8) 



Rx) 



1 

2x 



1 



1 ( ■sl2xf{x) + \+a- Vc^TTu V^ 
Kx \ J2xf{x) + 1 +a+ Va2 -^ 1 J 



-1 



(2.9) 
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(K^O, Ci and yo are constants). 



Eq. ( |2.6| ) provides a good example. The symmetry group is large enough to reduce to quadra- 
tures. This however really means that we have replaced a differential equation for y(x) by a func- 
tional equation (|2.9[) for f(x). To obtain a graph y = y(x) we still have to do numerical calculations. 



2.2 The difference invariants 

Let us consider four points xji^ on a line and the values yt = yix^) at these points : 

(x„_i , X„, X„+i , X„+2> 3'n-l > yn,yn+l, yn+l) 



(2.10) 



The SL(2,R) group generated by the prolongation of the vector fields p.l[ ) to the points p.lO| ) 
will transform these points in the (x,y) plane but will leave certain functions of them invariant. We 
calculate these invariants using known methods [2,. ..,7]. The result is that out of these coordinates 
we can construct five 4-point difference invariants 



jn _ yn yn-\ jn+l _ yn+l yn 

yXnXn—\ 



'\^n+ 1 ^n 



yXfi+lXn-l 



j„+2 _ yn+2 - yn+l 

1 ~ / '■ 

y Xn+2Xn+l 



jn+2 _ yn+2 - yn 

'2 ~ I 

y^n+2^n 



(2.11) 
(2.12) 



We mention that /""^\ /""^^ are just upshifts of /", I'^'^^ is an upshift of I!^'^^. Moreover, /", 7""^' 
and /j"*"^ involve coordinates of the first three points only. 



Let us now obtain difference schemes for the ODEs ( |2.4[ ) and ( |2.6| ). We put 

"« ~ ^n ~ ^n-lt "n+1 ~ ^n+l ~ ^nt '^n+2 ~ ^n+2 ~ ^n+l 

and expand yn+k = y{xn+k) about some point xq. First of all, we notice that the equation 
7^+1 -l'l = 



provides a good lattice. Indeed expanding (2.14i about the point x„ we obtain 



7«+l _ /« = y_^hn^^ _ h,) + ^^L_^^hi^ + hi) + ... 

In the continuous limit we put 

hn = a;„€ 



(2.13) 



(2.14) 



(2.15) 



(2.16) 



where a„ are constants of the order a„ ~ 1 and take e — > 0. From eq. ( |2.14[ ) and ( |2.15[ ), we see 
that for e — > we have 



On+i -an- 0(e) 
(if ^^ is finite) 



(2.17) 
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Let us now approximate the continuous invariants /j and I^ of eq. ( |2.3[ ). To obtain /^ , we need 
3 points (n - \,n,n+ \) ox {n,n + l,n + 2). From the expansion of 7""*^ and Z""*"^, we see that the 
correct expression is 



rn+l 






2xy" + y' 



y 



,/3 



+ 2{hn+, - h„)x \ J^ + O(h^) (2.18) 



3y 



/4 



and similarly for J^^^ (an upshift of /;■+'). We can now consider II1,I1*\I1^^,J1*\J1+^] as a 
basis for the four point difference invariants. Both 7"^^ and 7"^^ have the differential invariant /j 



as their continuous limit. Moreover, on the lattice (2.14i, this is an approximation of the order e 



(see(2.15l, (2.171, (2.181) 



To approximate /^ we must use all 4 points. Indeed, we have 



T / jn+2 _ jn+l 
^n+2 _ -' I -^l -^l 



2[ii + ii+^ +i'l+^) y 



r rrr o //2 



= -^{y y - 3^ ) + (/j„+2 - hn)R(x,y,y ,y ,y ,/ 



' ii" i/" -M('^)^ 



(2.19) 



where R is some differential expression (that is easy to calculate). Thus K"^^ goes to I^ for 
6^0 and provides an approximation of order e^ on the lattice ( 2.14 ) (since ( 2.14| ) implies also 
jn+i ^ jn+i =in = i^ with h independent of n). 



2.3 Difference schemes for the second order equation 

We put 



rn+i 



r. 



W1+1 — f'l — J 



(2.20) 



From eq. ( |2.20[ ) we obtain 



IT' =h{\ll + 2)^13 . 



(2.21) 



Since I\ is a constant (depending on the initial conditions (xo,3'o, x\,y\)) /3 will also be a con- 
stant. 



The scheme ( |2.20| ) can be solved explicitly for Xn+\ and j„+i and we obtain 



Xf,+ \ — Xfi-\ 



yn - yn-i 



PXn-l -(yn-yn-\) 



yn+Y 



P^n-iyn - (yn - yn~\)yn-\ 
PXn-\ -{yn-yn-l) 



(2.22) 



Thus, the invariant scheme is an explicit and linear one and moreover it is of the order e^. By 
comparison, a standard scheme of order e^ will be implicit and y„+i will be obtained by solving 
a cubic equation. An explicit standard scheme will be of order e. We mention here that it is the 
difference scheme that converges to the ODE like e^. This does not guarantee that the same is true 
for the solutions. 
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2.4 Difference scheme for the third order ODE 



From eq. p. 19 1 we see that an invariant difference scheme for the third order equation ( 2^ i is 



obtained by putting 

^„+2 ^ /7(/j'+l), /«+2 = /«+l = jn = ; _ (-2.23) 

Alternatively, we can put 



™+l ^ jn+2 



f^n+2 ^ p I _\ _^ I ^ jn-,2 = jn = J _ (2.24) 



Both schemes converge to the ODE like e^. The scheme ( |2.23[ ) can be solved explicitly for 

Xn+2, yn+2 and we obtain 

^n yn ~ (^n+iyn+1 ,~ ~cr\ 

Xn+2 ^ -^ yn+2 ^ ; (2.25) 

where we have defined 

1 n 7 17/- rn+l\ 1 r«+ln 



a;„+i--^[2/if(7«^0 + /r]AH^ • (2-26) 

4 V ^n+l 

We stress that I\ is a constant (independent of n and determined by the initial conditions) but 
7"^' and hence w„+i depend on n and must be calculated at each step using 
4 

^r' ^ —^=^—^\yn+\ - yn-\ - 2/l ^JXn+lXn-l\ (2.27) 

Thus the SL(2,]R) invariant scheme is a very simple one for any function F{Ji): x„+2 and yn+2 
are obtained explicitly in terms of their values at n - 1 , n and n + I. The standard scheme will 
be nonlinear and hence implicit. The condition /" = /""*"' in ( |2.23| ) is actually quite restrictive for 



third (and higher) order ODEs. The constant I is determined by the initial conditions and p.23| ) 
imposes a relation between y(0), y'iO) and ^"(0). For general initial conditions a relation of the 
type /""""' = yl" is more suitable (y is determined in terms of the initial conditions). 

3 One dimensional realization of sl(2,R) 

3.1 The Lie algebra and differential invariants 

The Lie algebra is realized by the vector fields 

Xi=dy, X2=ydy, X^=y^dy (3.1) 

so X is an invariant quatity. The lowest order differential invariant is the Schwartzian derivative 

h = ^^(y'y"'-\y"^) ■ (3.2) 

Thus the equation 

^2^y'y"' - ^/") = Fix) (3.3) 



will be invariant under this realization of SL(2,]R) for any function F{x). For F{x) ^const. eq ( |3.3| ) 
is invariant under GL(2,]R) generated by ( |3.1[ ) and X4 - dx- For F{x) - the ODE is invariant 
under SL(2,R)(g)SL(2,R). 
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3.2 Invariant difference scheme 



The ODE ( |3.3[ ) can be approximated by a four point difference scheme. The space of difference 
invariants is five dimensional and is generated by 



r,n+2 _ (yn+2 " yn)CV»+l " Jn-l) 

t^ — , , , , ! ^n-\t ^nt ^n+\t ^n+2 



{yn+2 - yn+i)(yn+\ -yn) 
and the ODE p.3|) is approximated by 



A _ pn+2 

hn+i ^ h„ = h, F{xn,0) ^ F{x) 



(3.4) 

(3.5) 
(3.6) 



We have chosen a uniform lattice but any equation E{Xfi-\,Xn, x„+i, x„+2) = will provide an 
alternative invariant lattice. 



The scheme p.5| ) can be solved explicitly and we obtain 

Kniyn - yn-\)yn+\ " (jn+l " yn-\)yn 



Xn+2 = Xn+\+h = {n + 2)h + Xq , yn+2 

where we have put 

/j2 / h 
Kn=4 l-yFU + - 



K„{y„ -yn-i)- {y„+i -y„-i) 



Expanding y""*"^ in (|3.5|) about the point jc„ + ^ we obtain 



Kn-^,{y'y"'-^y"') + 0{h^) 



(3.7) 



(3.8) 



(3.9) 



Thus the scheme ( |3.7[ ) approximates the ODE ( |3.3[ ) with second order accuracy and the scheme 
is explicit and linear. 

A standard scheme which is second order accurate will be implicit and yn+2 will be calculated 
from an algebraic equation (at least quadratic). 



4 Numerical analysis 

In this section we concentrate on the two-dimensional realization of sl(2,l.) and apply the symme- 



try preserving and standard schemes to ( |2.4[ ) and ( |2.5[ ). For some numerical results concerning the 
one-dimensional realization of sl(2,R) see Ref.[ 1]. 



4.1 Second order equations 

The general solution of ( |2.4[ ) is given in ( |2.7[ ) for C i^ 0. On Fig|T]we show the exact solutions yi 
(increasing branch) and y2 (decreasing branch) for y = 150, yb = 5,C = e^. The step for the exact 
solution was h = 0.05. The symmetry preserving method integrates yi up to the singularity at 
X = xq = C/y ~ 20.3 then continues along the second branch j2 to its 'initial' value. The standard 
method fails to converge close to the singularity (where the solution becomes complex). 



Symmetry preserving discretization of SL(2,]R.) invariant equations 9 

4.2 Third order equation 

Let us now consider eq. (|2.6|) with a = -I. We put I"'''^ = yl" where y is determined by the initial 



'1 /'I 

conditions. As the step h tends to zero, y will tend to y = 1 as required in eq. ( |2.14| ). 



On Fig|2]we compare the accuracies of the standard and symmetry preserving scheme. Since 
in this case no exact analytic solution is available we compare with a reference solution using a 
Matlab Runge-Kutta scheme with a tolerance on the error set at tol= 10"^. The initial conditions 
were set at yo = 1» Jq - 10, y'^ = -4 and this corresponds to a solution with no singularity on the 
real axis for < x < 16. We see that the accuracy is better for the symmetry preserving scheme 
by a factor of 10. 

A singular solution is shown on Figj3j The initial conditions were set at jo = 1> Jq ~ ^' ^o ~ ^ 
and a singularity occurs for x ~ 1.7. Matlab solvers and standard schemes stop providing solutions 
close to the singularity. The symmetry preserving method approaches the singularity closely and 
then continues along the second branch of the singular solution towards an appropriate initial 



condition. Qualitatively we have the same features as for the second order equation ( |2.4| ). The 
solution itself stays finite but its derivative becomes infinite at the singularity. 



5 Conclusions 

The main conclusion that we draw from the examples so far considered is that the symmetry pre- 
serving schemes are actually simpler than the standard ones. We have shown that they can be 
explicitly solved (see eq. ( |2.22[ ) and ( 2.25 1) while preserving second order accuracy. Not sur- 



prisingly, the numerical calculations confirm that the symmetry preserving schemes have better 
accuracy and provide better results close to singularities. 

A specific conclusion from Fig {T] and Fig|3} born out by other examples [1], is that the sym- 
metry preserving method provides qualitative information about singularities of the solutions, not 
available from standard methods. We can not only pinpoint the position of the singularity. We also 
see from the curves that the solution itself remains finite, but its derivative becomes infinite. In 
addition we see that it is a square root type singularity. Indeed if we follow the solution backwards 



from the singularity at x = Xh we see two different branches. For the second order equation (2.4 1 



this behaviour is obvious from the explicit solution ( |2.7| ). For the third order equation ( |2.6| ) this is 
not visible without numerical calculations. 
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Figure 1. Behaviour of the symmetry preserving scheme near the singularity for eq. (|2.4| 
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Figure 2. Discretization errors for standard and symmetry preserving schemes for eq. ( |2.6| l, a - -I for a 
regular solution 
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Figure 3. Behaviour of the symmetry preserving scheme near a singularity for eq. (|2.6| 



